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^ Abstract 

Numerical solvers of the incompressible Navier-Stokes equations have reproduced turbulence phenomena 

such as the law of the wall, the dependence of turbulence intensities on the Reynolds number, and 
experimentally observed properties of turbulence energy production. In this article, we begin a sequence 
of investigations whose eventual aim is to derive and implement numerical solvers that can reach higher 

^ Reynolds numbers than is currently possible. Every time step of a Navier-Stokes solver in effect solves 

JY^ ^ a linear boundary value problem. The use of Green's functions leads to numerical solvers which are 

^ highly accurate in resolving the boundary layer, which is a source of delicate but exceedingly important 

physical effects at high Reynolds numbers. The use of Green's functions brings with it a need for careful 
quadrature rules and a reconsideration of time steppers. We derive and implement Green's function 
based solvers for the channel flow and plane Couette flow geometries. The solvers are validated by 
reproducing turbulent signals which are in good qualitative and quantitative agreement with experiment. 

^ 1 Introduction 

T—i The incompressible Navier-Stokes equations are given by du/dt + (u.V)u = —Vp + Au/Re, 

^ where u is the velocity field and p is pressure. The incompressibility constraint is V.u = 0. We 

assume that a characteristic speed U and a characteristic length L have been chosen and that 
the Reynolds number Re is given by UL/u, where i' is the kinematic viscosity. It is assumed 

CN that the unit for mass is chosen so that the fluid has a constant density equal to 1. 

^ The topic of this paper is the use of Green's functions to solve the incompressible Navier- 

^ Stokes equations. The Navier-Stokes equations arc nonlinear while Green's functions are based 

on linear superposition. Thus the solutions of the incompressible Navier-Stokes equations 
^ cannot be described using Green's functions. However, if we discretize the Navier-Stokes 

equations in time but not in space, and the time discretization treats the nonlinear advection 
term (u.V)u explicitly and the pressure term — Vp and the viscous diffusion term Au/Re 
implicitly, each time step is a linear boundary value problem. The simplest such discretization, 
which is to treat the advection term using forward Euler and the pressure and diffusion term 
using backward Euler, gives the equation 

. + ((u.V)u)" = -Vp"+^ + -^Au"+i 
At Re 
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where the superscripts indicate the time step. This is a hnear boundary value problem for 
u""*"^ with the constraint V-u""*""^ = and with boundary conditions on u depending upon the 
geometry of the flow. Green's functions may be derived for this linear boundary value problem 
as shown in the theory of hydrodynamic potentials |16] . 

Green's functions exploit the principle of linear superposition to express solutions of linear 
boundary value problems in integral form. In numerical methods based on Green's functions, 
the weight of the method falls upon quadrature rules as opposed to rules for the discretization of 
derivatives. The importance of quadrature rules is already clear in the early work of Rokhlin 
\2H\ I29| . where Richardson extrapolation and trapezoidal rules are used to effect accurate 
quadrature of the integral equations of acoustic scattering and potential theory. 

From the beginning, Greengard, Rokhlin, and others jlUI E51 ES] have emphasized the 
ability of Green's function based methods to handle very thin boundary layers. Shear flows 
such as channel flow or pipe flow or plane Couette flow are characterized by very thin boundary 
layers at high Reynolds numbers. There is turbulence activity in the boundary layer as well as 
in the outer flow and the viscous effects propagate into the domain from the boundary layer. 
Green's function based methods are likely to be advantageous to handle such boundary layers. 
It is legitimate to ask why a numerical method must be believed to capture the effect of the 
viscous term with the very small 1/Re coefficient. In Green's function based methods, that 
effect is captured exactly by the analytic form of the Green's function. 

As far as we are aware, time integration using Green's functions has not been tried on a 
nonlinear problem of the scale and difficulty associated with fully developed turbulence. Thus 
some of the issues that come up in relation to time integration in Section 3 cannot be considered 
unexpected. 

Many of the subtleties associated with the numerical integration of the Navier-Stokes equa- 
tions are related to the treatment of pressure. The equations do not explicitly determine the 
evolution of pressure. Instead, pressure is determined implicitly through the incompressibil- 
ity constraint on the velocity field. One of the key algorithms for solving the Navier-Stokes 
equations in channel and plane Couette geometries is due to Kleiser and Schumann 1980]. 
Kleiser and Schumann introduced a numerical technique for enforcing the physically correct 
boundary conditions on pressure. Another method was introduced by Kim, Moin, and Moser 
\14:\ 1987] in a paper that is a landmark in the modern development of fluid mechanics. Kim et 
al. reproduced several features of fully developed turbulence from direct numerical simulation 
of the Navier-Stokes equations. Their calculation was initialized with a velocity field that was 
generated using large eddy simulation. One of the highlights of the paper by Kim et al. is the 
correction of a calibration error in a published experiment using numerical data. 

The channel geometry is rectangular with x, y, and z being the streamwise, wall-normal, 
and spanwise directions by convention. The corresponding components of the velocity field u 
are denoted u, v, and w. The walls are at y = ±1 with fiuid in between. The no-slip boundary 
conditions require u = Oaty = ibl. The boundary conditions in the wall-parallel directions 
are typically periodic in numerical work. The fiow is driven either by maintaining a constant 
mass fiux or a constant pressure gradient in the streamwise direction. In plane Couette fiow, 
the geometry is the same but the walls are moving. The no-slip boundary conditions are 
{u, V, w) = (0, ±1,0) at y = ±1. Plane Couette fiow is driven by the motion of the walls. 

Both the Kleiser-Schumann and Kim-Moin-Moser methods come down to solving linear 
boundary value problems in the y or wall-normal directions. The periodic directions are tackled 
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Fig. 1.1: Very thin boundary layer at y = 1 in the solution of the fourth order boundary value 
problem (D^ — f3'^){D'^ — a'^)u = 1 with boundary conditions u{±l) = u'{±l) = and 
with parameters a = 10^ and /3 = ^/2a. The solid markers are from an exact formula 
and the solid line is a numerical solution. 



using Fourier analysis and dealiasing of the nonlinear advection term. Each Fourier component 
then yields a linear boundary value problem in the y direction. The y direction is discretized 
using Chebyshev points yj = cos j'tt/M with j = 0, . . . , M. 

Here we parenthetically mention the interpretation of the parameters a and j3, which 
occur in the ensuing discussion. The parameters are given by = ^^/A^ + n^/A^ and = 

+ J Re/ At, where i and n are the Fourier modes and 27: Ax and 2ttAz are the dimensions of 
the domain in the streamwise and spanwise directions, respectively. The parameter 7 depends 
on the time integration scheme. More details are found in Section 3. 

In Figure we have shown the solution of the linear boundary value problem 

{D' - f){D^ - a^)u[y) = f{y) n(±l) = n'(±l) = (1.1) 

with f{y) = 1. Here D = A fourth order boundary value problem of this type occurs 
explicitly in the method of Kim-Moin-Moser but it is treated as a composition of two second 
order boundary value problems corresponding to the factors — and — In the 
method of Kleiser-Schumann, a fourth order boundary value problem is not formed explicitly. 
Both methods solve the second order boundary value problem 

{D''-f)u{y) = !{y) u{±\) = ^ (1.2) 

by using the Chebyshev series n(y) = J2m=o '^rnTmiu) and the set of equations obeyed by the 
coefficients Cn given on page 119 of Gottlieb and Orszag ^8, 1977]. 

Although the method on p. 119 of Gottlieb and Orszag [8] has been extensively used in 
turbulence computations for more than two decades, its numerical properties have not been 
investigated as far as we know. For reliable use in solving the Navier-Stokes equations at high 
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Reynolds numbers, the method should be able to accurately reproduce thin boundary layers, 
such as the one shown in Figure [TTTj There is reason to be concerned. If the method is used to 
solve fourth order problems of the type (1.1), it forms linear systems with condition numbers 
of the order for the fourth order boundary value problem (1.1) and of order o? for the 

second order boundary value problem (1.2) [35j . For the problem shown in Figure [lT] the 
condition number is more than 10^^ and greatly exceeds the machine epsilon of double precision 
arithmetic. 

Zebib |ini 1984] and Greengard [S] suggested using a Chebyshev series for the highest 
derivative. For the fourth order problem (1.1), the Chebyshev expansion would be vi!'" = 
J2m=o '^rn.Tmiy) ■ This device avoids the ill-condition in Chebyshev differentiation due to clus- 
tering at the end points that causes large errors in spectral differentiation. However, the 
spectral integration method of Zebib and Greengard also has a condition number greater than 
for the fourth order boundary value problem (1.1 ) [H [35] . 

The method of spectral integration has been extended and investigated carefully by Viswanath 
|35j . He finds that the equations presented somewhat tersely on p. 119 of Gottlieb and Orszag 
[S] is in fact a form of spectral integration. The methods used by Kleiser-Schumann and 
Kim-Moin-Moser have numerical properties that are practically identical to that of Zebib and 
Greengard. The method of Zebib and Greengard has no advantages. This fact has not al- 
ways been appreciated. For example, Lundbladh, Henningson, and others [THl [T^] and Waleffe 
j36| implemented the Zebib-Greengard form of spectral integration believing it to be advan- 
tageous. When we refer to spectral integration, it includes the methods of Gottlieb-Orszag, 
Kleiser-Schumann, Zebib, Kim-Moin-Moser, and Greengard as well as the more general and 
powerful formulation derived in ^35j. 

Regardless of which version of spectral integration is used, the fact remains that the linear 
system for the fourth order boundary value problem (1.1) has a condition number of . 
Yet, remarkably, even systems with condition numbers exceeding 10^^ (see Figure 1.1) can be 
solved with a loss of only five or six digits of accuracy. The accuracy of spectral integration in 
spite of large condition numbers can be partly explained using the singular value decomposition 
j35| . Another property of spectral integration (in all its forms) is that some of the intermediate 
quantities have large errors which cancel in the final answer [35]. A robust implementation 
must take these two properties into account. Spectral integration can indeed handle thin 
boundary layers, such as the one shown in Figure |1.1[ in spite of large condition numbers. 
The robustness of spectral integration was essential to the outstanding success of the methods 
of Kleiser-Schumann and Kim-Moin-Moser in more than two decades of use (however, not all 
implementations are equal). 

In Figure [lT] the thickness of the boundary layer is of the order 10"^. It takes more than 
10,000 Chebyshev points in the interval — 1 < y < 1 to resolve that boundary layer in spite 
of quadratic clustering near the endpoints. That is a lot more than the number needed if the 
grid points are chosen in a suitably adaptive manner. Viswanath [35 has derived a version 
of spectral integration that applies to piecewise Chebyshev grid points. Using that method, 
the number of grid points needed to solve a linear boundary value with a boundary layer as 
thin as the one shown in Figure 1.1 is reduced from 8192 to 96. It appears that this new 
method can be used to obtain considerable improvement in both the Kleiser-Schumann and 
Kim-Moin-Moser methods. 

The Green's function method, whose development we begin in this paper, is an alternative 
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which in its final form will enjoy the same advantages. Spectral integration is an essentially 
one dimensional idea and cannot be generalized to pipe flows with non-circular cross-sections 
and to other non-rectangular geometries. The use of Green's functions on the other hand will 
generalize. A great many analytic and numerical complications arise when Green's functions 
are derived for cross-sections of pipes as a part of a numerical method for solving the Navier- 
Stokes equations. It is essential to develop the method for the channel geometry, as we do here, 
before those difficulties are confronted. It is also possible that the Green's function method 
will turn out to be faster than the methods of Kleiser-Schumann and Kim-Moin-Moser, revised 
in the manner suggested in the previous paragraph, but one cannot be certain until the two 
alternatives are developed to their final form. Spectral integration over piecewise Chebyshev 
grids appears to be sensitive to the location of the nodes used to divided the interval [35] . The 
Green's function method is likely to be much less sensitive. Lastly, we mention that the Green's 
function method has a theoretical advantage. When implemented using suitable quadrature 
rules, its numerical stability is immediately obvious. 

The Green's function method for solving the Navier-Stokes equations in channel and plane 
Couette fiow geometries is developed in Sections 2 and 3. The quadrature rule that is used 
is provisional. The way to derive robust quadrature rules is indicated in Section 3 and the 
complete method will be given in the sequel to this paper. 

In Section 4, we show the theoretically intriguing result that Green's functions may be 
used to eliminate numerical differentiation in the wall-normal or y direction entirely from 
the numerical scheme. The derivatives that occur in the nonlinear advection term can be 
transferred to the Green's function using integration by parts with the result that numerical 
differentiation is replaced by analytic differentiation. Such a scheme is not practical at high 
Reynolds numbers for reasons given in that section. 

In Section 5, we validate the Green's function based method. In view of the extensions 
discussed in this introduction, the code has been written in such a way that it can reach 
hundreds of millions of grid points using only a dozen or two processor cores. Since the 
piecewise Chebyshev extension with robust quadrature rules is yet to be fully developed, the 
full capabilities of this code are not exercised. Yet we report simulations with up to ten million 
grid points and investigate certain aspects of fully developed turbulence to demonstrate the 
viability of the Green's function approach. 



2 Green's functions and template boundary value problems 

Every time step in the solution of the Navier-Stokes equations in the channel geometry reduces 



to the solution of a number of linear boundary value problems of the type (1.1) and (1.2). In 
this section, we derive the Green's functions for the solutions of those boundary value problems. 
The Green's functions can be derived using very standard methods. However, the resulting 
expressions are unsuitable for numerical evaluation. When the parameters a and (3 are as 
large as 10^, as in Figure [lT] quantities of the type e^* or e~^* will overflow. Thus we begin 
by deriving the Green's functions in a manner that leads to expressions suitable for accurate 
numerical evaluation. In the last part of this section, we consider the evaluation of derivatives 



such as du/dy, where u is the solution of either of the boundary value problems ( 1.1 ) and (|1.2 



and the evaluation of the solution u when the source term / is given in the form / = dfi/dy. 
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2.1 Green's functions of linear boundary value problems 

Let Lu = M^"^ + ai(y)u^"'~^'* + • • • + a„_i(y)n'^"'^^ + a„(y)ti. The coefficients ai(y), 1 < « < n, are 
assumed to be real-valued and sufficiently smooth. The adjoint operator is given by L^v = 
(— l)"^^"^ + (— l)"~"^(aiu)*^"~"'^^ + • • • + anV. We assume throughout that the functions that 
arise have the requisite order of smoothness and that n>2. The degree of differentiabihty is 
specifically mentioned only if there is a nontrivial reason for doing so. 

The lemmas in this subsection are not new. They can be found in [5] in some form or the 
other. However, our derivation leads to formulas which are easier to manipulate and which 
are suitable for numerical evaluation. Our derivation of the Green's function for the boundary 
value problem Lu = f, a < y < c, with suitable boundary conditions on u, is based on the 
Lagrange identity, which is the next lemma. 



Lemma 1. For any two functions u and v, the Lagrange identity v Lu — uL~^v 
with 

n— 1 n— fc— 1 
k=0 r=0 



[uv]' holds, 



□ 



and ao = 1. 

Proof. Begin with J v Lu dy and integrate each term by parts repeatedly. 
Define 

/ u \ 

uW 



The quantity [uv] which appears in the Lagrange identity may be written as [uv] = uT'Av, 
where A is an n x n matrix. All the entries of A are determined by the lemma. However, all 
that we need to know about A is that it has the following reverse triangular structure 



A 



\l 



(-1)"-! \ 



and that the reverse diagonal is as shown above. 

Let ui,. . . ,Un be a basis of solutions of the homogeneous problem Lu = 0. Similarly, let 
vi, . . . ,Vn be a basis of solutions of the adjoint problem L~^u = 0. Denote that n x n matrices 

and (wi, . . . , {;„) 

by U and V, respectively (the determinant of U is the Wronskian). The Lagrange identity 
implies the following lemma. 



Lemma 2. 



dy 



U^AV] = 0. 
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We will assume that the bases of solutions are chosen in such a way that 



U{yYA{y)V{y) = I 



{2.1) 



where / is the identity matrix. The homogeneous solutions Ui and Vi are used to construct the 
Green's function of Lu = f. Before deriving the Green's function, we give the identity 



/ 



det 



Ml 
U2 



/ 



det 



/ ui u'l 

U2 U'2 



(n-l) \ ■ 



(2.2) 



,{n-l) 

(n-l) 



The entry equal to 1 in the last column of the numerator is in row number j. This identity is 



derived as follows. We choose the j-th column of (2.1 ) to get 

/ \ 



/ 



U{y)^A{y) 



\ 



(n-l) 



1 



V / 



Because of the reverse triangular structure of A, the last entry of Avj is equal to Vj. Identity 



(2.2) is implied by Cramer's rule. By working with rows of (2.1) instead of columns, we get 



the identity 



det 



i-iy 



Vl 
V2 



,{n-2) qX 



(n-2) 



{n-2) 
Vn 



/ 



det 



Vl 
V2 

\ Vn 



Xn-1) 



(n-l) , 
Vn I 



(2.3) 



The identities (2.2) and (2.3) are used to construct Green's functions in Section 2.2. 



So far, we have not specified the boundary conditions u must satisfy in addition to Lu = /. 
We take the domain to be a < 7/ < c and require that u{(i) must lie in an n — £ dimensional 
subspace Vn (this corresponds to I linear conditions on n(a)). Similarly, the right boundary 
conditions require that u{c) should lie in a n — r dimensional subspace Yr- We require £+r = n. 
The Green's function is built up using homogeneous solutions of Lu = which satisfy the left 
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or the right boundary conditions. We assume that the basis solutions are chosen and then 
ordered in such a way that 



ui{c), . . . ,Ui{c) and ue+i{a) , . . . , Un{a) 



span the subspaces Vr and Vi, respectively. The following lemma gives the boundary conditions 
satisfied by vi, . . . ,Vn, a basis of solutions of L^v = which is related to ui, . . . , ii„ by (2.1 1. 



The lemma is useful for checking correctness of the implementation. It may also be used for 
the construction of Vi given Uj. Its proof is obvious from U'^AV = I. 

Lemma 3. vi{a), . . . ,vi{a) span the orthogonal complement of the n — i dimensional space 
A{a)'^Vi and ?);+i(c), . . . , Vn{c) span the orthogonal complement of the n — r dimensional space 
A{cfVr. 

Let u be the solution of Lu = f subject to the boundary conditions u{a) £ Vi and u{b) G Vr- 
If we apply the Lagrange identity using u and Vi, where Vi is a solution of the homogeneous 
problem L^v = 0, we get fvi = ^ (ji^Avij for i = 1, . . . , n. The equations with i = !,...,£ 
are integrated from a to y and the rest are integrated from y to c. The boundary conditions 
as given by the previous lemma imply that 



fvi 



u{y)'^ A{y)i)i{y) 



fvi 
fvi+i 



u{y)'^ A{y)vi{y) 
u{y)'^ A{y)vi+i{y) 



fVr, 



u{y)'^ A{y)vn{y). 



The last entry of A{y) u is equal to (— l)"'n. Using Cramer's rule, we get 



det 



(n-2) 



(n-2) 



det 



V2 v'2 



,in-l) X 



,("-!) 



(n-1) 
Vn 



The following lemma gives the Green's function in a more useful form. 



(2.4) 



Lemma 4. The solution of Lu = f subject to the boundary conditions u{a) G Vn and u{c) G Vr 
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is given by 

u{y) = ui{y) I vi{t)f{t)dt + --- + ui{y) I vi{t)f{t)dt 

J a 

rc re 
-ui+iiy) / Vl+i{t)f{t)dt un{y) / Vn{t)f{t)dt. (2.5) 



Proof. Use (2.3) and (2.4). 



□ 



The following lemma justifies the delta-function interpretation of Green's functions favored 
by physicists. It is used in Section 2.3 and Section 4. 

Lemma 5. If ui, . . . ,Un and vi, . . . ,Vn are bases of solutions of the homogeneous problems 
Lu = and L^v = 0, respectively, that are related by U (y)"^ A{y)V (y) = I, we have 



i=l 



and 



E 

i=l 



UiV. 



for i = 0,l,...,n-2 

1 for j = n — 1 

for j = 0,l,...,n-2 

(-!)"-! for j = n-l. 



Proof Use (2.2) and (2.3). □ 



2.2 Template boundary value problems 

The first template boundary value problem is 

[d^ -P^)u = f 

with boundary conditions ^(±1) = 0. The Green's function of this boundary value can be 
deduced in any number of ways. We have 

_ /3(-2+s/+t) , -/3(4~s,+t) I p/3(-y+i) _ -P(2+y+t) 

0(v. «) = ^ 2,i(e-» 1) 

for —l<t<y<l. The Green's function is symmetric and we have u{y) = G{y, t) f{t) dt. 
This form of the Green's function suits numerical work because none of the terms will overflow 
for even (3 very large. The terms in the numerator are factored as follows: 

g/3(-2+j/+t) ^ 

g-/3(4-j/+t) ^ g-2/3g-/3(l-y)g-/3(l+t) 
g-/3(2+j/+t) ^ ^-/3{l+y)g-/3(l+f)_ 
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None of the factors will overflow even for large /3. The term e^^"^"*"*^ is not factored and will not 
overflow because t < y. Using these factorizations and noting that y and t must be exchanged 
to get the Green's function for y < t, we infer that the evaluation of the solution of the first 
template problem using u = f^^ G{y, t) f{t) dt reduces to the quadratures 



-fiit+i) 



fit) dt, 



\-t'(t+^)f{t)dt, r e''^-^+''>f{t)dt, [\''^-^+'^f{t)dt (2.7) 



and 



-1 



"^J'\y-t\ 



y 



fit) dt 



(2.8) 



with fj, = 13. Each one of these quadratures yields a function of y and must be multiplied 
by a prefactor which is also a function y. For example, the first term of (2.6) contributes the 
_g/3(-i+j/)/2/3(e~'^^ — 1) to e^^~^^^\f(t) dt. Since the term is unchanged when y 



prefactor —e^^~^~^y^ /2l3{e~^^ — 1) to f^^e^^~^^^^ f{t) dt. Since the term is unchanged when 
and t are interchanged to obtain the Green's function in y < t region, it contributes the same 
prefactor to e^^~^^^^ f{t) dt. The prefactor of the function defined by (2.8) with 7 = /3 is 

l/2^(e~^'' — 1). Unlike the result of the quadratures (2.7) and (2.8), the prefactors do not 



depend upon / and can be computed and stored in advance. Thus the cost of solving the 



first template problem is very nearly equal to the cost of the quadratures (2.7) and (2.8) with 
/3. 
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The second template problem is the fourth order boundary value problem 



j u 



f 



with boundary conditions u{±l) = ti'(ibl) = 0. For this problem, it takes more work to write 
the Green's function in such a way that there are no numerical overflows even if a and /3 are 
very large. However, the final result is similar to what we have seen for the first template 
problem. The evaluation of u can be reduced to the quadratures (2.7) and (2.8) with 7 = a 
and ^ = p. The results of quadratures are multiplied by prefactors and summed to obtain u. 
A certain 4x4 matrix we will now derive is useful for computing the prefactors. 

Like the first template problem, the second template problem is self-adjoint. We take the 
basis of homogeneous solutions to be 



Ul 






(y-i) 


_ g-"(j/-i) 




U2 




ae-^(^-i)- 


-/3e° 




-a(y-l) 


Us 






(y+i) 






U4 


= ae'^(^+i)- 


ae-^(^+^)- 


-/3e" 


(^+i)+/3e- 


-a{y+l) 



It may be verified that ui and U2 satisfy the right boundary conditions while M3 and U4 satisfy 
the left boundary conditions as assumed in Section 2.1. The functions vi,V2,V3,V4 may be 
calculated using (2.2) or Lemma |3j It is convenient to define 



W 



2_2„-4a 



-4a/3(5VV 



4q/35W^^ + 4 



4a 



+ 32 a^(3'^6^e 



a 135^ 
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a S^a^ 



13 a 



-2 ^-2 c 



13 a 



-4 ;9 a (5 (7 " + 4 /3 Q, 5 o- e-2 /3 
-2/3(52 + 2/3 52e-4/5 
-8/32a^e-2/3 -2 '^^'"""/"'° +2/3^a2e-2" 



13 a 5^ 



Pad-' 



8/32a<5e-2'3-2« + 2 ^ - 2/35a2e-4^ 
-2/352 +2^52e-'*'' 



Fig. 2.1: Entries of a 4 x 4 matrix with the first two columns placed above the last two. 

This matrix determines the Green's function of the second template problem. Here 
6 = - 13"^ and a = a + {3. 

where 5 = — jS"^ and cr = q + /3 (the Wronskian is equal to e^'^'^'^^W). The function vi is 
equal to 

2 a /3 (/3 + a) (-/3 + af e" (^"^^ - 2 a /3 (/3 + a) {-P + af ^^''^'^ +2a (3 {-(3 + a) {P + af e"" (^+3^ 
+2al3 i-P + a) (/3 + a)^ e"'^ (^+3)+2 a /3 (/3 + a) (-/3 + a)^ e-'*"-''^-3^-2 a /3 + a) {p + af 6-^"+^ 
-4a2/3 (-/? + a) (/3 + q) e-^"-^ ?'-/5_2 a/3 (/3 + a) + af e-^^-°'y"^"'-2 a^ {-P + a) (/3 + af e-4/3+"2/~" 
-4 a 13'^ (-/3 + a) (/3 + a) e'^ z^"" ^'-"+4 a /^^ (-/3 + a) (/3 + a) e" ^"^ ^^+4 a^^j (_;3 + (/3 + a) e'^ ""^ ''"^ « 

divided by W . The expression for V2 is similarly long. 

For t < y, the Green's function is given by G{y,t) = ui{y)vi{t) + U2{y)v2{t). This Green's 
function is determined by the 4x4 matrix shown in Figure |2.1[ We think of the rows and 
columns of the matrix as being indexed by — /?, /3, —a, a in that order. The (— /3, — /3) entry, 
which appears in the top left corner, is divided by W to get the coefficient of e~^(^"'"^-'e~'^^*"'"^'' 
in the expression for G{y, t) for t < y. The other entries are interpreted similarly but there 
are two special entries. These are the (— entry which must be interpreted as W times 
the coefficient of e^^*"^-* and the (—a, a) entry which must be interpreted as W times the 
coefficient of e"^*"^-*. The Green's function for y < t is obtained from symmetry. 

It follows that solving the second template problem (D^ — (D^ — a^) u = f with bound- 



ary conditions u{±l) = u'{±l) = reduces to quadratures (|2.7|) and (2.8) with 7 = and 



7 = /3. Each quadrature yields a function of y which is multiplied by a prefactor. The prefactor 



is determined using the 4x4 matrix of Figure 2.1 and the formula for W. 



The formula for W and the entries of the 4x4 matrix use 5 and a to avoid cancellation 
errors. Because of the way the parameters a and /? arise in the numerical integration of channel 
flow, 5 = — can be evaluated accurately. 



Each of the quadratures (2.7) and (2.8) is well- conditioned. However, if a ~ /3 there will 
be large cancellation errors when the results of the quadratures are multiplied by prefactors 
and summed. This phenomenon may be understood as follows. When a ^ /3, the solutions 
gitaj/^ g±/3j/ fQj.jjj basis of homogeneous solutions. When a = /3, the basis is e^"^, ye^°^. When 
a ~ /3, the Green's function tries to produce terms which resemble ye°^ using terms such as 
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(e"y - ef^y^ /{a- (3) resulting in large cancellation errors. Fortunately, this situation does not 
arise in channel flow or plane Couette flow. 

2.3 Derivatives using Green's functions 

For the template boundary value problems Lu = f, we have derived Green's functions such 
that u{y) = J\G{y,t) f{t)dt. Here we will consider the use of Green's functions to evaluate 
derivatives such as du/dy. 

The ability to differentiate solutions of boundary value problems using Green's functions 
has been utilized in an important paper by Greengard and Rokhlin [10^. They consider the 
boundary value problem u" + p{y)u' + q{y)u = f{y) and solve it by representing the solution u 
in the form u = f^^ G{y, t) a{t) dt, where G is the Green's function of a linear boundary value 
problem with constant coefficients which satisfies the same boundary conditions. In fact, the 
background boundary value problem is simply taken to be u" = f. With the representation of 
u using a, the boundary value problem becomes an integral equation for a. Starr and Rokhlin 
|31| have generalized the method to first order systems. The papers by Greengard, Rokhlin, 
and Starr show how to apply numerical methods based on Green's functions to problems with 
non-constant coefficients. Once the boundary value problem is cast in integral form using the 
background Green's function, the method handles diagonal blocks using Nystrom integration 
and pieces together the global solution efficiently by exploiting the low rank of the off-diagonal 
blocks. 

We derive integral formulas for derivatives of solutions of both the second and fourth order 
boundary value problems. In addition, we consider boundary value problems of the type 
Lu = dfi/dy and Lu = and show how to get the solution u as well as its derivatives 

without numerically differentiating /i or /2. In Section 4, these calculations are used to show 
that numerical differentiation in the wall-normal or y direction can be entirely eliminated in 
the numerical integration of channel flow. 

For the template second order problem, which is Lu = d'^u/dy'^ — f3'^u = f with boundaries 
tt(±l) = 0, we take the Green's function to be G{y,t) = ui{y)vi{t) for — 1 < t < y < 1. The 
Green's function for — 1 < y < i < 1 is taken to be G{t,y) since the problem is symmetric 
or self-adjoint. The function ui satisfies the right boundary condition and the relationship 
between ui, M2 and ^1,^2 is as given in Section 2.1. The Green's function for — l<y<t<lis 
also given by —U2{y)v2{t). As a consequence of symmetry, we have —U2{t)v2{y) = ui{y)vi{t). 

If G{y,t) = ui{y)vi{t), we have Gi{y,y) = G2{y,y) + 1, where subscripts of G denote 
partials with respect to the first or second argument. This follows from symmetry and Lemma 
|5j In addition, we have G{l,y) = G{y,—1) = because ui satisfies the right boundary 
condition and vi satisfies the left boundary condition. 

The solution of Lu = / is given by 

u{y) = G{y, t) fit) dt + 1^ G{t, y) f{t) dt. (2.9) 
Differentiating with respect to y, we get 

u'iy) = J"^ Gi{y, t) fit) dt + J^' G2{t, y) f{t) dt (2.10) 
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where subscripts of G stand for partial differentiation. The integral equation is no longer 
symmetric in y and t. Suppose the boundary value problem is Lu = dfi/dy. We may substitute 
/{ for / in (2.9) and integrate by parts to get 

u{y) = - J''^G2{y,t) flit) dt- J^' Gi{t,y) flit) dt. (2.11) 

This integral equation for uiy) is not symmetric. Differentiating with respect to y and using 
Giiy, y) = G2iy, y) + 1, we get 



dv ry 

-^ = - j Gi2iy,t)fiit)dt- j Gi2it,y)fiit)dt + fiiy). (2.12) 



The template fourth order problem is Lu = (D^ — /3^)(L'^ — a'^)u = f with boundary 
conditions ui±l) = u'i±l) = 0. We again take the Green's function to be Giy,t) for —1 < 
t < y < I- From Section 2.1, we have G(y,t) = ni(?/)fi(t) + U2iy)v2it). Figure 2.1 gives 



the coefficients of the Green's function as explained in Section 2.2. The functions ui(y) and 
U2iy) satisfy the right boundary condition. Using symmetry, we take the Green's function for 
— l<y<t<ltohe Git, y). As a consequence of symmetry, we have 

uiiy)viit) + U2iy)v2it) = -U3it)v3iy) - U4(t)f4(y)- 

Using this identity and Lemma |5] we deduce that 

Giiy,y) = G2iy,y) 
Giiiy,y) = G22iy,y) 
Giiiiy,y) = G222iy,y) + 1- 

Here the subscripts of G denote partial differentiation with 1 and 2 standing for the first and 
second arguments of G. Since ui and U2 satisfy the right boundary conditions while vi and V2 
satisfy the left boundary conditions, we have 

G(l,y) = Giil,y) = G(y, -1) = G2(y, -1) = 0. 

In other words, the Green's function satisfies the boundary conditions. 

The solution of Lu = /, with the boundary conditions associated with the fourth order 
template problem, are given by 

uiy) = r Giy, t) fit) dt + C Git, y) fit) dt (2.13) 



J-l Jy 

as before. Differentiating with respect to y gives 

dy 

d^i fy ^ 

dy'^ 



l'^Giiy,t)fit)dt + 1^ G2it,y)fit)dt (2.14) 
j'Giiiy,t)fit)dt + J' G22it,y)fit)dt. (2.15) 
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The subscripts of G denote differentiation. If the fourth order template problem is of the form 
Lu = dfi/dy, its solution is given by 



G2(y,t)/i(t) dt 



Giit,y)h{t)dt. 



(2.16) 



This form of the solution is obtained after substituting dfi/dt for / in (2.13) and then inte- 
grating by parts. The boundary terms vanish. Differentiating with respect to y, we get 

du '•^ 



dy 
d^ 
dy'^ 



Gi2{y,t)fi{t)dt- 
Gn2{y,t)fi{t) dt 



Gi2{t,y)h{t)dt. 



Gi22{t,y)h{t) dt. 



(2.17) 
(2.18) 



The boundary terms vanish on both occasions. If the template fourth order boundary value 
problem is in the form Lu = d'^f2/dy'^, the analogous formulas are as follows: 



u{y) 

du 
dy 
<Pu 
dy^ 



G22{y,t)f2{t)dt+ / Gu{t,y)f2{t)dt 



-1 
y 



Gl22{y,t)f2{t)dt+ / Gu2it,y)f2it)dt 



Gll22{y,t)f2{t)dt + 



-1 



G'ii22(t, y) f2{t) dt + (Gi22(y, y) - G^iy, y)) f2{y). 

(2.19) 

These formulas are derived using the properties of G given in the previous paragraph. 

Formulas (2.9) through (2.19) give a method to compute solutions and solution derivatives 
without numerical differentiation even when the right hand side of the boundary value problem 
is given as a derivative. If the right hand is dfi/dy, these formulas use /i and not dfi/dy. The 
derivatives are transferred to the Green's function which can be differentiated analytically. In 
the case of the template fourth order problem, the kernels of the formulas can be described 
using a matrix such as the one displayed in Figure 2.1 In fact the kernels can be obtained 
by multiplying the entries of that matrix with suitable powers of a and /3. With such a rep- 
resentation the kernels can be evaluated in a numerically stable way as described in Section 
2.2. 

The numerical evaluation of formulas (2.9) through (2.19) is affected by discretization and 
rounding errors in varying ways. To avoid writing down long formulas, we limit the discussion 
of numerical errors to the template second order problem and note that very similar issues 
arise for the template fourth order boundary value problem. 

When the integral formulation is used, the solution of the template second order problem 
(-D^ — /3^)ti = / at the boundary point y = 1 is obtained as the sum of the following four terms: 



1 e^(*-i)/(t) 



2^(1 



-4/3^ 



dt, 



e2^^V^(*+i)/(t) 
_i 2/3(1 -e-4/3) 



1 g-2/3g-/3(t+l) ^(^) 



2/3(1 



1 p/3{t-l) 



m 



2/3(1 



dt, I __4^^ dt, I ^ ^/^(^ _^-4i3'^dt. 

(2.20) 

The second and third terms are exceedingly small even for moderate /3. The main contribution 
to numerical error is from the exact cancellation between the first and the last terms. The 



2 Green's functions and template boundary value problems 



15 



magnitude of the first or the last term is of the order \f\aQ / f^"^- If the quadrature rule is a very 
good one, each of the integrals may be evaluated with an error of around \f\^P~'^emachine- 
If such a quadrature rule is devised, the error in the boundary layer will also be of the same 
order. If we suppose / = /3^, then the exact formula will look like 1 — e^^^~^^ near the y = 1. 
Since 1 is a special number in machine arithmetic, the subtraction y — 1 will be exact at y = 1 
but not at other nearby points. If we take the subtraction error at y = 1 to follow the same 
model as at other points, we get the error in the boundary layer using the exact formula to 



be of the order 



1 



/3e or l/l^ f-machine/ So we See that the use the integral form has 



the potential to be more accurate in the boundary layer than even the mathematically exact 



formula. Here we envisage quadrature rules for the sort of integrals that occur in (2.20) which 



take into account the occurrence of terms such as e ^(*+^) in the integrands and whose weights 
and nodes are computed using extended precision. 



The use of formulas such as (2.10) to compute the derivative du/dy is especially accurate 



in the boundary layer. For instance, at y = 1 the first term of (2.20) gets multiplied by j3 and 



the last term by —(5 with the result that there is no cancellation error in the boundary layer. 
In view of this observation, some of the errors reported in Table 7 of [TU] may appear a little 
high for the function derivative. 



Finally, we consider a type of numerical error that arises in formulas such as (2.16) that 
express the solution of Lu = dfi/dy in integral form and without differentiation of /i. If /3 is 
large in the template second order problem (D^ — I3'^)u = /, the solution satisfies u ~ —f/P"^ 
away from the boundary. Thus if / is given in the form dfi/dy, the solution will satisfy 



u 



dy I 



-(5 and a formula such as (2.16) essentially has to produce the derivative of /i away 



from the boundary using integration. Differentiation is defined by subtracting nearly equal 
quantities and the cancellation errors inherent in that process cannot go away entirely. The 



same comment applies to formulas such as (2.10) which produce solution derivatives using 



an integral formula or to the evaluation of solution derivatives using the background Green's 
function as in [TU] or to the method of spectral integration discussed in the introduction. The 
principle contribution to the solution of (D^ — /3^)u = / for large /3 and away from the boundary 
is due to the term 

1 



p\^-y\ fit) dt. 



1 



2/3(e-4/3 - 1) 

This is the term which makes the solution look like — away from the boundary. When 
(2.10) is used to evaluate du/dy with u being the solution of Lu = /, the leading contribution 



is from the two terms 



2^(e-4/3 - 1) 



fit) dt 



.I3{y-t) 



fit) dt 



Here the cancellation error we are looking for is evident. A numerical method that uses integral 
formulas to evaluate solution derivatives or to eliminate derivatives that appear on the right 
hand side would benefit by treating such terms together, especially when quadrature rules are 
derived. 
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3 Time integration of the Navier-Stokes equations 

Let u = (u, w) be the velocity field of channel flow or plane Couette flow. We assume the 
domain to be periodic in the wall-parallel directions with periods equal to 2ttAx and IttAz 
in the streamwise and spanwise directions, respectively. The Fourier decomposition of the 
velocity field is given by 

L/2 N/2 

E E u,,„(y)e^^-/^^+-^/^^ 

e=-L/2n=N/2 

This Fourier decomposition assumes the number grid points in the streamwise and spanwise 
directions to be L and A^. The notation u^_„ denotes a Fourier coefficient of the entire velocity 
field. Similarly, Uf^„, V£^ni we,n denote the Fourier coefficients of the streamwise, wall-normal, 
and spanwise components of the velocity field, respectively. The components of the vorticity 
V X u are denoted by ^, 77, C and their Fourier components are denoted similarly. 

Often which modes i and n apply is clear from context and the Fourier modes are indicated 
as u, V, w without subscripts. The i = n = modes are the mean modes and are denoted 
using an over-bar. For example, the mean mode of the streamwise velocity is u. In both the 
flows considered here, the range of the y variable is — 1 < y < 1, with the walls located at 
y = ±l. 



3.1 The Kim-Moin-Moser equations 

We take the Navier-Stokes equations to be du/dt+H = —Vp+Au/Re, with H = {Hi, H2, H^) 
being the nonlinear term. Both the Kleiser- Schumann [T5]and Kim-Moin-Moser methods 
begin by substituting the truncated Fourier expansion of the velocity field u. The various 
Fourier modes are coupled through the nonlinear term. The nonlinear term is dealiased using 
the 3/2 rule 

Both the methods use identical equations for the mean streamwise velocity and mean 
spanwise velocity: 



du - 1 d'^u dw - 1 d'^w 

- = -Hi+p, + —^ — = -H, + —^ (3.1) 

For plane Couette flow pg = and the boundary conditions are u(±l) = ±1 and id{±l) = 0. 
For channel flow, u(±l) = tD(ibl) = but pg is non-zero. We may take pg = 2/ Re and 
maintain a constant pressure gradient or we may take 



1 /■+^ . 1 du 



1 - 
P9 = 2j^ ^^^y 



2Re dy 



(3.2) 



and keep the streamwise mass flux ^ X^l^ u dy constant at 2/3. The laminar solution of channel 
flow is u = (1 — y^, 0, 0) in both cases. 
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n? 




A2 




V 




e 


n? 




A| 



The equations for the (£, n) mode are 

Here ah the hatted variables are Fourier coefficients of the (£, n) mode and are functions of y. 
As before D denotes d/dy. The incompressibihty constraint V.u = gives i£u/Ax + dv/dy + 
inw/ h.z = 0. The equations are solved in this form by the Kleiser-Schumann method. In the 
Kim-Moin-Moser method these equations are altered to 

'Al"Al 




i- 














Re 


AS 


"Af 


± 




e 


r? 


Re 


AS 


~ A2 



^^^^ = irAD-^-T-. \ ^- (3-3) 



The boundary conditions are ?}(±1) = ^{±1) = ^(±1) = 0. Here 

dz dx ^ dx^ dz^ dydx dydz 



The entire velocity field can be recovered using u, w, r], and v [1] 

Imposing physically correct boundary conditions on pressure causes some complications 
and is a potential pitfall. Early discussions of this issue are found in |15| |21] . A thorough 
discussion of this topic, important both for mathematical theory and for computation, is found 
in an illuminating paper by Rempfer [27 \ . 

3.2 Time stepping using Green's functions 



If the Kim-Moin-Moser equations (3.1 ) and (3.3) are discretized in time, we get linear boundary 



value problems in the wall-normal or y direction. Green's functions will be used to solve these 
linear boundary value problems. An advantage of this method is that the boundary layers are 
analytically built into the Green's functions. 

The original paper by Kim, Moin, and Moser pij used the CNAB (Crank-Nicolson and 



Adam-Bashforth) discretization in time. If the method is applied to the fj equation in (3.3), 
we get 

At 2 ^R~e[ "Al"A|jl^ 2 j" 

The superscripts denote time steps. It is well-known that the numerical stability of Crank- 
Nicolson can be dicey in spite of its stability region being the entire left half plane. The 
eigenvalue equal to A corresponds to an amplification factor of (1 -|- AAt)/(l — AAt). The 
amplification is by a factor less than 1 in magnitude for eigenvalues with a negative real part. 
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However, the amplification factor can be very close to 1 for eigenvalues such as A = —10^° 
which correspond to rapid decay. If care is taken to use the same scheme for differentiating 
f/"^"^ and 77", or if the boundary value problem is solved for f^^+i -i-f^" at each time step, CNAB 
be stable. We found it difficult to stabilize CNAB for the v equation in (3.3). This could be 
because we are mixing integration using a Green's function with the second derivative that 
comes from the left hand side of (3.3), or it could be because the best possible quadrature 
rules for this problem are yet to be derived. We will not consider CNAB any further. 

Suppose dX/dt = f{X)+AX/Re, where f{X) is a nonlinear term. The time discretizations 
we consider are of the following form: 



1 

At 



j=0 



j=0 



If the u equation of (3.1) is discretized in time, it fits the template second order boundary 
value problem (D^ — I3'^)u = f with 



u = u 



n+l 



At 



(3.4) 



The time discretization of the w equation of (3.1) fits the template second order boundary 
value problem (D^ — /3'^)m = / with 



~Af' 



s-l 



s-1 



(3.5) 



The time discretization of the fj equation of (3.3) also fits the template second order boundary 
value problem: 



^2 

^ Al^ Ay At 



Re 



s-l 



(3.6) 



i=o 



The time discretization of the v equation of (3.3) fits the template fourth order boundary value 
problem (1.1) as follows: 



-,n+l 



a 



+ 



n 



AV 



2 , 

" + AT' 



/ 



Re 



At ^ 



a ]v 



+ ReJ2b,H; 



(3.7) 

From the manner in which and f3^ arise, the advantage of casting the Green's function for 
the template fourth order problem using 5 = — /J^ and a = + /J^ , as we did in Figure 



2.1 and Section 2.2, is evident. Both 6 and a can be evaluated without cancellation errors. 



The equations (3.4), (3.5), (3.6), and (3.7) are reduced to quadratures of the type (2.7) and 
(2.8) as explained in Section 2. In the fourth order problem (3.7), which is solved for t)""^^, 
the right hand side uses second derivatives of v"'~^ from previous stages. The calculation of 
this second derivative is reduced to quadratures using (2.14). 
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Fig. 3.1: Absolute stability regions of the explicit halves of the implicit-explicit methods based 
on backward Euler, CNAB, BDF2, and BDF3. 



The time stepping schemes we have implemented use 

s = 1, 7 = 1, ai = -1, 5i = 1 
s = 2, 7 = 3/2, ai = -2, as = 1/2, bi = 2, 62 = -1 
s = 2, 7 = 11/6, ai = -3, a2 = 3/2, 03 = -1/3, 61 = 3, 62 = -3, 63 = 1- 

These are implicit-explicit multistep schemes that correspond to backward Euler, BDF2, and 
BDF3 respectively. For derivations of these schemes, see jU [3 

The absolute stability regions of the explicit halves of some time integration schemes are 



shown in Figure 3.1 In turbulence simulations, the nonlinear advection term, which is dis- 
cretized using an explicit scheme, is more of a constraint on the time step than the diffusion 
term which is handled implicitly. The discretization of the the viscous diffusion term is by itself 
unconditionally stable. As explained at length by Rempfer [27^, the reason for treating the 
diffusion and pressure terms implicitly is to get the physically correct boundary conditions on 
pressure. Thus the constraint on the time step is determined mainly by the absolute stability 
regions of the explicit part. Since the nonlinear term corresponds to advection, absolute sta- 
bility regions which run along the imaginary axis should be preferred. In this respect. Figure 



3.1 shows that BDF3 is superior to the other methods. 

To complete the description of these methods, we need to explain the method is used 



to solve the quadrature problems (2.7) and (2.8) numerically. These quadrature problems 



are extremely well-conditioned even for large fi. The method that is currently implemented 



for (2.7) is to expand the integrand in a Chebyshev series and integrate the terms of the 
series using well-known formulas. A better method would be to obtain quadrature nodes and 
weights for weighted integrals with weight functions equal to e^^^*^^-*, evaluate the other factor 
f(t) at the nodes using an accurate and efficient interpolation algorithm, and sum using the 
quadrature weights. Such a method will be developed in future research. The current method 
develops spurious difficulties when fi is large, although it is good enough to allow us to exhibit 
simulations of fully developed turbulence in Section 5. In addition, if the idea of representing 
functions using piecewise Chebyshev collocation, which is briefly mentioned in the introduction 
and discussed at greater length in the context of spectral integration in [35j , is employed even 
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the basic quadrature that is now implemented is hkely to be adequate, even for very large 
The Green's functions of Section 2 are completely independent of the discretization used in the 
wall-normal or y direction. The discretization could be Chebyshev, or piecewise Chebyshev, or 
something else. The ease with which piecewise Chebyshev discretization can be incorporated 
into numerical methods that use Green's functions was one of our prime motivations. At the 



moment, the quadrature problem (2.8) is solved numerically using spectral integration. Similar 



comments apply to this quadrature problem as well. 

3.3 Why is time integration of the Navier-Stokes equations stable? 

The nonlinear term in the Navier-Stokes equations is usually given in advection form: H = 
(u.V)u. It can be recast into the rotation form H = a; x u or the divergence form H = V.(uu) 
or the skew-symmetric form H = V.(uu) /2 -|- (u.V)u/2. The form in which the nonlinear term 
is implemented numerically is known to influence the stability of the time stepping scheme. 
An early discussion based on the centered-difFerence discretizations is given in [20] . When the 
nonlinear term is calculated with dealiasing, as we do here, it is usual to employ the rotation 
form. However, when the nonlinear terms are not dealiased, the skew symmetric form appears 
to have some advantages [illl |39] . 

We cannot expect an easy answer to why the time integration of Navier-Stokes is numeri- 
cally stable. Numerical stability theory deals with simple linear equations such as du = \u and 
ut = Uxx- Even though the stability theories are an indispensable guide to practical computa- 
tion, it must be kept in mind that they are based on equations whose mathematical stability is 
trivial. Since there is no mathematical stability theory of the 3D Navier-Stokes equations, we 
cannot expect to understand the numerical stability of its discretizations perfectly. Practical 
computation shows that stability can be quite subtle. During one of our simulations using the 
method of Kim, Moin and Moser l^Aj and BDF3 time discretization, the computation became 
unstable after more than a million time steps. 

Here we suggest that energy cascades from low modes to high modes may be crucial to the 
numerical stability of time integration. The idea of energy cascades appears to have been first 
introduced by Obukhov [25j. It is often mistakenly attributed to Kolmogorov along with "the 
k~^^^ law." For a historically correct discussion, see [37j. The idea is that in a fully developed 
turbulent flow energy is injected at the lowest spatial modes and flows down to the finest scales 
where it is dissipated by viscosity. Viscosity is deemed to be unimportant in the intermediate 
inertial range. Shear flows are neither homogenous nor isotropic, especially near the walls. Yet 
in both channel flow and plane Couette flow, energy is injected to the mean mode via the u 



equation in (3.1). A glance at the equations for v and fj (3.3) shows that the equations are 
dominated by dissipation if £ and n are large. For such fine scales, numerical stability will not 
be an issue. Thus it appears possible that numerical time integration is stable because energy 
flows from the mean mode to the finest scales, where viscous effects dominate the nonlinear 
term and dissipate energy. As evidence in support of this hypothesis, we mention the familiar 
fact that under-resolved simulations show greater turbulence and a simulation which is severely 
under-resolved will become unstable. 
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4 A discrete model without spatial difFerentiation 



Here we explain how numerical differentiation in the wall-normal or y direction can be com- 
pletely eliminated by employing the divergence form of the nonlinear term. 



The Kim-Moin-Moser equation for u given in (3.1) has an Hi term. The nonlinear term 

dxiu"^) + dy{uv) + dz{uw 



So we may take Hi 



Hi is given by Hi 

u is advanced in time by solving the boundary value problem (3.4). The right hand side of the 



dyuv. The mean mode 



template boundary value problem is taken to be /, where / is given by (3.4). The Hi terms in 



that right hand side may be removed and a new right hand side written as dfi/dy introduced 
in their place. The contribution of a given time step to fi is taken to be Re times uv evaluated 
at that time step. These contributions are weighted by bj and combined as before. The mean 
mode u is advanced in time using (2.9) and (2.16) and its derivative, if needed in (3.2), is 



calculated using (2.10) and (2.12). 



The w and fj equations are treated similarly. The fj equation in (3.3) has an H4 term. The 
terms of H4 that do not require differentiation in y are 



+ 



92 



uw 



uw 



2„„2 



d^w 



dxdz dz^ 



dx^ dxdz 



and the terms which require a single differentiation are 

d'^uv d'^vw 
dydz dxdy 



These terms are separated and some of them are removed from the / given in (3.6) and a new 
term dfi/dy is inserted in the right hand side. 

The treatment of -0 modes is a bit more elaborate. The right hand side H^ may be written 



as 



d^{uv) ^ d^{uv) ^ d^{vw) ^ d^{vw 
dx^ 



+ 



dxdz"^ ~^ dx'^d 



+ 



dz^ 



I dx'^dy 



dydz"^ 




uw] 




dxdydz dydz^ 

d^{vw) \ 
dy'^dz j ' 



where terms are grouped depending upon whether they require zero, one, or two differentiations 
with respect to y. The right hand side of the template fourth order problem which is given as 
/ in (3.7) may be rewritten as / -|- dfi/dy + d'^f2/dy'^, with none of /, fi, and /2 involving 



differentiation with respect to y. The integral equations ( 2.13| ) through (2.19) may be used 
to produce v, dv/dy, and even d'^v/dy'^. The first derivative dv/dy is needed when the full 
velocity field is reconstructed from u, w, and the modes fj and v [14 . Turning this into a 
practical method hinges on numerical issues discussed at the end of Section 2. Eliminating 
numerical differentiation with respect to y may be useful if large number of Chebyshev points 
are used in the y direction. However, it appears that piecewise Chebyshev grids can resolve 
boundary layers and internal while using a small number of Chebyshev points in each sub- 
interval. Handling piecewise Chebyshev grids after reducing each time step to quadratures of 



the form (2.7) and (2.8) is as each as f^ = + f^- piecewise Chebyshev grids, numerical 
errors due to differentiation are not a cause for concern. 
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Fig. 5.1: The first plot shows the variation of the mean shear at the wall as a function of time. 

The second plot shows the evolution of the mean flow during transition to turbulence. 
Units are given in terms of channel half-width h and laminar center-line velocity U. 



5 Numerical validation 

The numerical method described in Section 3 for solving the incompressible Navier-Stokes 
equations has been implemented and tested in numerous ways. Many earlier computations of 
plane Couette flow and channel flow have been reproduced with precision. In this section, we 
describe a few computations of fully developed turbulence. All the computations described 
here are for channel flow. Channel flow is used far more often than plane Couette flow in 
turbulence simulations. 

A useful summary of turbulence computations using channel flow is given by Toh and 
Itano [32] . The Reynolds number Re by itself is not a good metric to assess the difficulty of a 
turbulence computation because simple solutions such as the laminar solution can be computed 
easily at any Reynolds number. The metric must take into account both the Reynolds number 
and the kind of solutions that the simulation generates. One useful metric is obtained by 
taking the time average of du/dy at the walls, where u is the mean streamwise velocity, and 
then computing Rcr = \/Re x \du/dy\. The frictional Reynolds number Rcr is a good measure 
of the difficulty of the simulation. The highest Re^ reached appears to be 2000 in the work of 
Hoyas and Jimenez [1 2 . The lowest Rer at which one still observe turbulence appears to be 
around 100 [13j. Nikitin [231 121] has derived a method for solving the incompressible Navier- 
Stokes equations in orthogonal curvilinear coordinates. Nikitin's method uses staggered grids, 
centered differences, cell averages for nonlinear terms, and explicit projections to enforce the 
incompressibility condition. The same program can handle channel, pipe, eccentric pipe and 
other geometries. Nikitin's method has been used to simulate fully developed turbulence at 
Re^ of 500. 



Figure 5.1 shows a simulation of channel flow with Re = 10"^, A^; = 2.0 and A^ = 1.0. The 



grid parameters used L = 64, M = 128, and N = 64. The boundary condition used was pg 



2 /Re in the equation for u given in (3.1 ). It is noticeable that the mean shear converges to —2 
at the upper wall. The equation for the mean flow u is given by du/dt + Hi = Pg + ■^d'^u/dy'^. 
The mean flow fluctuates very little once the flow is fully turbulent. If we average the mean 
flow over time, we get d'^u/dy'^ = Re{Hi — Pg) with ^^(±1) = 0. The Green's function for this 
boundary value problem is given by G{y,t) = {y — l)(t + l)/2 for — 1 < t < y < 1. Using 
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Fig. 5.2: Similar to Figure 5.1 but using a boundary condition that fixes the mass flux. 



(2.10), we get 



du 
di 



y=+l 



' Gi {y, t) fit) dt = ReC dt = -2. 

-l TT 



Since = 2/ Re, we have f_-^Hi{y)dy = 0. The reason Hi satisfies this condition appears 
not to be known. In rotating channel flows, the mean flow exhibits a stretch where its slope 
is given by the rate of rotation (see Figure 3 of Yang and Wu [38j). The reason for that 
phenomenon too appears to be unknown. The mean flow flattens during transition as evident 



from the second plot of Figure 5.1 At the very beginning, the mean flow develops oscillations 



which look somewhat like the oscillatory shears considered in fT7]. 



Figure 5.2 shows a turbulence simulation which flxed the mass-flux using (3.2 ). The param- 
eters used were Re = 10^, A^; = 1.0 and A^ = 0.5. Similar parameters are used by Moser et al. 
|22| . The grid parameters used L = 256, M = 256, and N = 128 correspond to approximately 
8.5 million grid points. This simulation reaches an Rer which is slightly less than 400. The 
mean flow again shows a tendency to oscillate during transition. 



When reporting further computations, the parameters used in Figure 5.2 will denoted as 
MKM. The other set of parameters we use is Re = 4000, A^^ = 2.0 and A^ = 1.0. Similar 
parameters were used by Kim, Moin and Moser [13] • Therefore this choice will be denoted 
KMM. The grid parameters L = M = N = 128 correspond to approximately two million grid 
points. 



Figure 5.3 shows the mean flow in the near- wall region for KMM and MKM. The flts to 
the law of the wall are quite good and of the same quality as in |14| [22] [52] . The variables 
with the superscript -|- are in frictional or wall units. 



Figure 5.4 shows the streamwise velocity recorded at various distances from the wall. Sig- 
nals recorded experimentally using hot wire anemometry look very much like the ones seen 



in the figure. For early experimental recordings, see [2]. Figure 5.4 is inspired more directly 
by the plots in [J. Hotwire anemometry remains the principal technique for obtaining quan- 
titative information about turbulent fiows, especially in the boundary layer. Even a quantity 
such as energy dissipation per unit mass, which is of much importance in turbulence theories, 
has to be obtained indirectly from turbulent signals of the type seen in the figure [SD] . Close 
observation of Figure 5.4 reveals that fluctuations at ~ 5, 10 are slightly more pronounced 
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Fig. 5.3: Fits to the law of the wall for KMM and MKM, respectively. The solid lines are 
from numerical simulation. The dotted lines in the viscous region are u"*" = in 
both plots. The dotted lines in logarithmic region are = 2.5 logy"*" + 5.5 and 
u'^ = 2.5 logy"'" + 5, respectively. 




tU/h X f tU/h 



Fig. 5.4: Turbulent signals for KMM and MKM, respectively. The vertical axis is u^{t) + 

and the signals correspond to y"*" ~ 5,10,20,30,40,50 in both plots. The factor / in 
KMM plot is taken to be 7.2/15 to make the time axes in the plots comparable. 
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for MKM than for KMM indicating the effect of the outer flow on near wah fluctuations. 
This phenomenon is somewhat surprising because the means fit the same curve in the viscous 
sublayer. It was recognized by interpreting experimental data recorded at various Reynolds' 
numbers . Examples of signals in simulations of isotropic turbulence are found in [33] . 

Figure |5.4| may lead us to suspect that the turbulent signals are not smooth. Such an 
appearance of non-smoothness is because of the fact that the fluctuations are governed by 
many different time scales. Figure 5.5 blows up a signal to reveal that the signal is in fact very 
smooth. 



6 Acknowledgments 

We thank Sergei Chernyshenko and Nikolay Nikitin for helpful discussions. This research was 
partially supported by NSF grants DMS-0715510, DMS-1115277, and SCREMS-1026317. 



References 

[1] U.M. Ascher, S.J. Ruuth, and B.T.R. Wetton. Implicit-explicit methods for time- 
dependent partial differential equations. SIAM Journal on Numerical Analysis, 32:797- 
823, 1995. 

[2] GK Batchelor and AA Townsend. The nature of turbulent motion at large wave-numbers. 
Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 
199(1057):238-255, 1949. 

[3] RF Blackwelder and RE Kaplan. On the wall structure of the turbulent boundary layer. 
Journal of Fluid Mechanics, 76(1):89-112, 1976. 

[4] J. P. Boyd. Chebyshev and Fourier spectral methods. Dover, 2001. 

[5] E.A. Coddington and N. Levinson. Theory of ordinary differential equations. McGraw- 
Hill, 1955. 



6 Acknowledgments 



26 



[6] E.A. Coutsias, T. Hagstrom, and D. Torres. An efficient spectral method for ordinary 
differential equations with rational function coefficients. Mathematics of Computation, 
65(214):611-636, 1996. 

[7] M. Crouzeix. Une methode multipas implicite-explicite pour I'approximation des equations 
d'evolution paraboliques. Numerische Mathematik, 35(3):257-276, 1980. 

[8] D. Gottlieb and S.A. Orszag. Numerical Analysis of Spectral Methods: Theory and Ap- 
plications. Society for Industrial and Applied Mathematics, 1993. 

[9] L. Greengard. Spectral integration and two-point boundary value problems. SI AM Journal 
on Numerical Analysis, 28:1071-1080, 1991. 

[10] L. Greengard and V. Rokhlin. On the numerical solution of two-point boundary value 
problems. Communications on Pure and Applied Mathematics, 44(4):419-452, 1991. 

[11] K. Horiuti. Comparison of conservative and rotational forms in large eddy simulation of 
turbulent channel flow. Journal of Computational Physics, 71(2):343-370, 1987. 

[12] S. Hoyas and J. Jimenez. Scaling of the velocity fluctuations in turbulent channels up to 
Rcr = 2003. Physics of fluids, 18:011702, 2006. 

[13] J. Jimenez and P. Moin. The minimal flow unit in near-wall turbulence. Journal of Fluid 
Mechanics, 225(213-240), 1991. 

[14] J. Kim, P. Moin, and R. Moser. Turbulence statistics in fully developed channel flow at 
low Reynolds number. Journal of Fluid Mechanics, 177(1):133-166, 1987. 

[15] L. Kleiser and U. Schumann. Treatment of incompressibility and boundary conditions in 

3-D numerical spectral simulations of plane channel flows. In 3rd Conference on Numerical 
Methods in Fluid Mechanics, volume 1, pages 165-173, 1980. 

[16] O. A. Ladyzhenskaya. The Mathematical Theory of Viscous Incompressible Flow, 2nd ed. 

Gordon and Breach, New York, 1969. 

[17] Y.C. Li and Z. Lin. A resolution of the Sommerfeld paradox. SIAM Journal on Mathe- 
matical Analysis, 43:1923, 2011. 

[18] A. Lundbladh, D.S. Henningson, and A.V. Johansson. An efficient spectral integration 
method for the solution of the Navier-Stokes equations. Technical Report FFA TN, 1992- 
28, 1992. 

[19] A. Lundbladh, D.S. Henningson, and S.C. Reddy. Threshold amplitudes for transition in 
channel flows. Transition, Turbulence and Combustion, pages 309-318, 1994. 

[20] N.N. Mansour, P. Moin, W.C. Reynolds, and J.H. Ferziger. Improved methods for large- 
eddy simulations of turbulence. In Symposium on Turbulent Shear Flows, volume 1, 
page 14, 1977. 

[21] P. Moin and J. Kim. On the numerical solution of time-dependent viscous incompressible 
fluid flows involving solid boundaries. Journal of Computational Physics, 35(3):381-392, 
1980. 



6 Acknowledgments 



27 



[22] R.D. Moscr, J. Kim, and N.N. Mansour. Direct numerical simulation of turbulent channel 
flow up to Rcr = 590. Physics of Fluids, 11:943, 1999. 

[23] N. Nikitin. Finite-difference method for incompressible Navier-Stokes equations in arbi- 
trary orthogonal curvilinear coordinates. Journal of Computational Physics, 217(2) :759- 
781, 2006. 

[24] N. Nikitin. Third-order-accurate semi-implicit Runge-Kutta scheme for incompress- 
ible Navier-Stokes equations. International Journal for Numerical Methods in Fluids, 
51(2):221-233, 2006. 

[25] AM Obukhov. Spectral energy distribution in a turbulent flow. Izv. Akad. Nauk SSSR, 
5(110):453, 1941. 

[26] K.N. Rao, R. Narasimha, and MA Narayanan. The bursting phenomenon in a turbulent 
boundary layer. Journal of Fluid Mechanics, 48(02):339-352, 1971. 

[27] D. Rempfer. On boundary conditions for incompressible Navier-Stokes problems. Applied 
Mechanics Reviews, 59:107, 2006. 

[28] V. Rokhlin. Solution of acoustic scattering problems by means of second kind integral 
equations. Wave Motion, 5(3):257-272, 1983. 

[29] V. Rokhlin. Rapid solution of integral equations of classical potential theory. Journal of 
Computational Physics, 60(2):187-207, 1985. 

[30] K.R. Sreenivasan and R.A. Antonia. The phenomenology of small-scale turbulence. Annual 
Review of Fluid Mechanics, 29(l):435-472, 1997. 

[31] P. Starr and V. Rokhlin. On the numerical solution of two-point boundary value problems 
II. Communications on Pure and Applied Mathematics, 47(8):1117-1159, 1994. 

[32] S. Toh and T. Itano. Interaction between a large-scale structure and near- wall structures 
in channel flow. Journal of Fluid Mechanics, 524:249-262, 2005. 

[33] L. Van Veen, S. Kida, and G. Kawahara. Periodic motion representing isotropic turbu- 
lence. Fluid Dynamics Research, 38(l):19-46, 2006. 

[34] J.M. Varah. Stability restrictions on second order, three level finite difference schemes for 
parabolic equations. SIAM Journal on Numerical Analysis, 17:300-309, 1980. 

[35] D. Viswanath. Spectral integration of linear boundary value problems, arxiv.org, 2012. 

[36] F. Waleffe. Homotopy of exact coherent structures in plane shear flows. Physics of Fluids, 
15:1517, 2003. 

[37] AM Yaglom. Alexander Mikhailovich Obukhov, 1918-1989. Boundary-Layer Meteorology, 
53(1), 1990. 

[38] Y.T. Yang and J.Z. Wu. Channel turbulence with spanwise rotation studied using helical 
wave decomposition. Journal of Fluid Mechanics, 692:137, 2012. 



6 Acknowledgments 



28 



[39] T.A. Zang. On the rotation and skew-symmetric forms for incompressible flow simulations. 

Applied Numerical Mathematics, 7(l):27-40, 1991. 

[40] A. Zebib. A Chebyshev method for the solution of boundary value problems. Journal of 
Computational Physics, 53(3):443-455, 1984. 



